c  This program reads the unformatted head file from the model output
c  and makes a list of the computed heads versus measured heads.
c  The measured heads must be in a file that contains layer,row,column
c  water-level elevation, and any description desired.  The number of
c  known water levels is currently limited to 200. The format of
c  the file must be (i5,i5,i5,f10.2,a90)
c
c  modified and converted to DG by Lynette Brooks
c  September 11, 1991
c
      PARAMETER LENX=40000
      real X(200,200),value,diff(400),total,abstotal,errsumsq
      integer arrow(400),arcol(400),arlay(400)
      real mean,absmean,resvar,dry
      integer count,maxcol,seqnum,rowcol
      CHARACTER FILENAME*96,ans*1,rest*90
      dimension text(4)
c  X array is the model head values
c  value is the known head value at a specified lay,row,col
c  diff is an array containing the difference between the modeled
c      value and the known value
c  arrow,arcol,arlay are arrays of the layer, row, and column
c  maxcol is the maximum number of columns in the model
c  total is total of differences
c  abstotal is total of absolute difference
c  errsumsq is error sum of squares
c  mean is mean of differences
c  absmean is mean of absolute differences
c  resvar is residual variance
c  dry is the value printed if the cell dries up during model run
c  count counts the number of known water levels
c  rest is description of known water levels
c  text(4) is needed to read unformatted file correctly
c
C OPEN NODE INTEREST FILE
   10 write(6,*)' ENTER NAME OF FILE CONTAINING KNOWN HEADS: '
      read (5,'(A)',ERR=10) FILENAME
      OPEN (unit=1,FILE=FILENAME,STATUS='OLD',ERR=10)
C OPEN UNFORMATTED INPUT FILE
   20 write(6,*)' ENTER NAME OF UNFORMATTED MODEL HEAD OUTPUT FILE: '
      read (5,'(A)',ERR=20) FILENAME
      OPEN (unit=2,FILE=FILENAME,STATUS='OLD',FORM='UNFORMATTED',
     CERR=20)
      REWIND (unit=2)
C OPEN OUTPUT FILE
   30 write(6,*)' ENTER NAME OF OUTPUT FILE: '
      read (5,'(A)',ERR=30) FILENAME
      OPEN (UNIT=7,FILE=FILENAME,STATUS='UNKNOWN',ERR=30)
c
c  prompt for value to use if node is dry
c
         write(6,*)'enter value for dry nodes'
         read(5,*)dry
c
c
C READ FIRST ARRAY RECORD FOR PROBLEM DIMENSIONS
c
      read (2) KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,KLAY
C ALLOCATE SPACE IN THE X array and stop program if not large enough
      NODES=NCOL*NROW
      write (6,1010) nodes,LENX
      IF (LENX.lt.nodes) then
         write (6,*)' Model array is too large.  Program halts'
         STOP
      end if   
c
c  rewind file so it can be read again
c
      REWIND (unit=2)
c
c  initialize the counter for number of comparisons
c
c  count = 0
c
c  does the user only want one time step,period,or layer
c  and if so, which one
c  set variables = 0 if user does not specify them
c
      istep = 0
      iper = 0
      ilay = 0
      write(6,*)' do you want only one time step, period, layer (y/n)?'
      read(5,'(a)')ans
      if (ans.eq.'y'.or.ans.eq.'Y') then
         write(6,*)' enter time step, period, layer'
         read(5,*)istep,iper,ilay
      end if
c   
c  READ A LAYER OF DATA
   50 read (2,END=70) KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,KLAY
      read (2,END=70) ((X(I,j),j=1,ncol),I=1,nrow)
      if (istep.ne.0) then     !istep ne 0 indicates user selected data wanted
         if ((kstp.ne.istep).or.(kper.ne.iper).or.(klay.ne.ilay)) then
            go to 50
         else
            go to 60
         end if
      end if
c
c  if user hasn't specified only one set of data, prompt for each one
c
      write(6,1020) kstp,kper,klay
      read(5,*)ans
      if (ans.eq.'n'.or.ans.eq.'N') go to 50
c  WRITE OUTPUT FILE FOR THIS LAYER if wanted
C REWIND NODE LIST FILE
      REWIND (unit=1)
c
c  print out header information
c
 60   write(7,1030)kper,kstp,klay
      write(7,*)' ROW  COL   KNOWN  CALCULATED  CALC-KNOWN  DESCRIPTION'
c
c  read known head values one line at a time, calculate difference
c  and write data
c
   62 read (1,'(3I5,f10.2,a90)',END=50) LAY,IROW,ICOL,VALUE,rest
      IF (LAY.NE.KLAY.OR.IROW.LT.1.OR.IROW.GT.NROW.OR.ICOL.LT.1.OR.
     CICOL.GT.NCOL) GO TO 62
      count = count + 1
      if (x(irow,icol).gt.1000000.or.x(irow,icol).lt.0.0001) then
         diff(count)= dry
      else
         diff(count) = x(irow,icol) - value
      end if
      arrow(count)=irow
      arcol(count)=icol
      arlay(count)=lay
      WRITE (7,'(2I5,3f10.2,a90)') IROW,ICOL,VALUE,
     1        x(irow,icol),diff(count),rest
      GO TO 62
c
c  all data read, do statistics
c  initialize totals
   70 total = 0.
      abstotal = 0.
      errsumsq = 0.
      write(6,*)' Number of Comparisons =         ',count
c
c  calculate totals
c
      do 80 i =1,count
         total = total + diff(i)
         abstotal = abstotal + abs(diff(i))
         errsumsq = errsumsq + (diff(i)**2)
 80   continue
c
c  calculate statistics and write results
c
      mean = total/count
      absmean = abstotal/count
      resvar = errsumsq/(count-2)
c
      write(7,*)' Number of Comparisons =         ',count
      write(7,*)' Total of Differences =          ',total
      write(7,*)' Total of Absolute Differences = ',abstotal
      write(7,*)' Mean of Differences =           ',mean
      write(7,*)' Mean of Absolute Differences =  ',absmean
      write(7,*)' Error Sum of Squares =          ',errsumsq
      write(7,*)' Residual Variance =             ',resvar
c
c  prompt for and open ascii file that can be read by info
c
      write(6,*)'do you want an ascii file for info?'
      read(5,'(a)')ans
      if (ans.eq.'y'.or.ans.eq.'Y') then

         write(6,*)'enter name of ascii info file'
         read(5,*)filename
         open (8,file=filename,status='unknown')
c
c  prompt for number of columns and do calculations and printout
c
         write(6,*)'enter number of columns in model'
         read(5,*)maxcol
         do 90 i=1,count
            seqnum = (arrow(i)-1)*maxcol + arcol(i)
            rowcol = 100*arrow(i)+arcol(i)
            idiff  = diff(i)
            write(8,1040)seqnum,rowcol,arlay(i),diff(i),idiff
 90      continue
      end if
 1010 FORMAT (I10,' ELEMENTS OF THE X ARRAY USED OUT OF',I10,
     -' POSSIBLE.')
 1020 format(' do you want data for time step',i5,' period',i5,
     1           ' and layer',i5,'(y/n)?')
 1030 format(' Data for Stress Period',i5,' Time Step',i5,
     1            ' Layer',i5)
 1040 format (i5,',',i5,',',i5,',',f10.2,','i10)
c
c  close files and stop
c
      CLOSE (unit=1)
      CLOSE (unit=2)
      ENDFILE (UNIT=7)
      CLOSE (UNIT=7)
      stop
      END
